library(data.table)
library(rmarkdown)
library(AnnotationHub)
library(tidyverse)
library(tximport)
library(ggplot2)
library(DESeq2)
library(pheatmap)
library(gridExtra)
library(ggplotify)
library(ggrepel)
library(UpSetR)
Assign your species of interest
AnnotationSpecies <- "Homo sapiens" # Assign your species
ah <- AnnotationHub(hub=getAnnotationHubOption("URL")) # Bring annotation DB
ahQuery <- query(ah, c("OrgDb", AnnotationSpecies)) # Filter annotation of interest
if (length(ahQuery) == 1) {
DBName <- names(ahQuery)
} else if (length(ahQuery) > 1) {
DBName <- names(ahQuery)[1]
} else {
print("You don't have a valid DB")
rmarkdown::render()
}
AnnoDb <- ah[[DBName]] # Store into an OrgDb object
# Explore your OrgDb object with following accessors:
# columns(AnnpDb)
# keytypes(AnnoDb)
# keys(AnnoDb, keytype=..)
# select(AnnoDb, keys=.., columns=.., keytype=...)
AnnoKey <- keys(AnnoDb, keytype="ENSEMBLTRANS")
# Note: Annotation has to be done with not genome but transcripts
AnnoDb <- select(AnnoDb,
AnnoKey,
keytype="ENSEMBLTRANS",
columns="SYMBOL")
# Check if your AnnoDb has been extracted and saved correctely
class(AnnoDb)
## [1] "data.frame"
head(AnnoDb)
## ENSEMBLTRANS SYMBOL
## 1 ENST00000543404 A2MP1
## 2 ENST00000566278 A2MP1
## 3 ENST00000545343 A2MP1
## 4 ENST00000544183 A2MP1
## 5 ENST00000286479 NAT2
## 6 ENST00000520116 NAT2
.sf files have been created from fastq data by salmon
# This code chunk needs to be written by yourself
# Define sample names
SampleNames <- c(paste0("A", 1:3), paste0("B", 1:3))
# Define group level
GroupLevel <- c("A", "B")
# Define contrast for DE analysis
Contrast <- c("Group", "B", "A")
# Define a vector for comparing TPM vs Counts effect
TvC <- c("TPM", "Counts")
# Define .sf file path
sf <- c(paste0("salmon_output/",
SampleNames,
".salmon_quant/quant.sf"))
# Define sample groups
group <- c(rep("A", 3), rep("B", 3))
# Create metadata
metadata <- data.frame(Sample=factor(SampleNames, levels=SampleNames),
Group=factor(group, levels=GroupLevel),
Path=sf)
rownames(metadata) <- SampleNames
# Explore the metadata
print(metadata)
## Sample Group Path
## A1 A1 A salmon_output/A1.salmon_quant/quant.sf
## A2 A2 A salmon_output/A2.salmon_quant/quant.sf
## A3 A3 A salmon_output/A3.salmon_quant/quant.sf
## B1 B1 B salmon_output/B1.salmon_quant/quant.sf
## B2 B2 B salmon_output/B2.salmon_quant/quant.sf
## B3 B3 B salmon_output/B3.salmon_quant/quant.sf
# Assign sample names to the input (.sf) file path
names(sf) <- SampleNames
# Run tximport
# tpm vs original counts
# input sf: a factor of all .sf files' path
txi_tpm <- tximport(sf,
type="salmon",
tx2gene=AnnoDb,
countsFromAbundance="lengthScaledTPM", # Extracts TPM
ignoreTxVersion=T)
txi_counts <- tximport(sf,
type="salmon",
tx2gene=AnnoDb,
ignoreTxVersion=T)
# tpm
head(txi_tpm$counts)
## A1 A2 A3 B1 B2 B3
## A2MP1 11.828289 9.640769 4.3214387 10.768983 9.930742 23.582072
## A3GALT2 0.000000 0.000000 0.9922462 4.934298 8.216794 9.121198
## A4GNT 3.878022 1.983684 0.9962060 5.005617 7.978542 4.061164
## AACSP1 99.547876 123.181812 54.2225647 15.620886 26.531886 22.692137
## AADACL2 0.000000 0.000000 0.0000000 1.016813 0.000000 0.000000
## AADACL4 1.951536 0.000000 0.9857014 0.000000 0.000000 0.000000
dim(txi_tpm$counts)
## [1] 12202 6
# counts
head(txi_counts$counts)
## A1 A2 A3 B1 B2 B3
## A2MP1 12.118 14.265 3 10.045 12.188 16.12
## A3GALT2 0.000 0.000 1 5.000 8.000 9.00
## A4GNT 4.000 2.000 1 5.000 8.000 4.00
## AACSP1 97.000 129.000 59 15.000 24.000 24.00
## AADACL2 0.000 0.000 0 1.000 0.000 0.00
## AADACL4 2.000 0.000 1 0.000 0.000 0.00
dim(txi_counts$counts)
## [1] 12202 6
# Set a function creating dds and vsd
dds_vsd_fn <- function(txi) {
# Create a DESeq object (so-calledd dds)
des <- DESeqDataSetFromTximport(txi,
colData=metadata,
design=~Group)
# Create a vsd object (so-called vsd)
ves <- vst(des, blind=T)
# Output them as a list
return(list(dds=des, vsd=ves))
}
TPM <- dds_vsd_fn(txi_tpm)
Counts <- dds_vsd_fn(txi_counts)
# Outputs
# dds: TPM/Counts[[1]] or TPM/Counts[['dds']]
# vsd: TPM/Counts[[2]] or TPM/Counts[['vsd']]
# TPM
TPM[[1]]
## class: DESeqDataSet
## dim: 12202 6
## metadata(1): version
## assays(1): counts
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(0):
## colnames(6): A1 A2 ... B2 B3
## colData names(3): Sample Group Path
head(counts(TPM[[1]]))
## A1 A2 A3 B1 B2 B3
## A2MP1 12 10 4 11 10 24
## A3GALT2 0 0 1 5 8 9
## A4GNT 4 2 1 5 8 4
## AACSP1 100 123 54 16 27 23
## AADACL2 0 0 0 1 0 0
## AADACL4 2 0 1 0 0 0
# Counts
Counts[[1]]
## class: DESeqDataSet
## dim: 12202 6
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(0):
## colnames(6): A1 A2 ... B2 B3
## colData names(3): Sample Group Path
head(counts(Counts[[1]]))
## A1 A2 A3 B1 B2 B3
## A2MP1 12 14 3 10 12 16
## A3GALT2 0 0 1 5 8 9
## A4GNT 4 2 1 5 8 4
## AACSP1 97 129 59 15 24 24
## AADACL2 0 0 0 1 0 0
## AADACL4 2 0 1 0 0 0
# TPM
TPM[[2]]
## class: DESeqTransform
## dim: 12202 6
## metadata(1): version
## assays(1): ''
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(4): baseMean baseVar allZero dispFit
## colnames(6): A1 A2 ... B2 B3
## colData names(4): Sample Group Path sizeFactor
# Counts
Counts[[2]]
## class: DESeqTransform
## dim: 12202 6
## metadata(1): version
## assays(1): ''
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(4): baseMean baseVar allZero dispFit
## colnames(6): A1 A2 ... B2 B3
## colData names(3): Sample Group Path
using ‘avgTxLength’ from assays(dds), correcting for library size gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates
# Set a function estimating size factors, dispersions, and perform wald test
DESeqPrep_fn <- function(List) {
List[[1]] <- estimateSizeFactors(List[[1]])
List[[1]] <- estimateDispersions(List[[1]])
List[[1]] <- nbinomWaldTest(List[[1]])
return(List)
}
# Update dds with the function
Counts <- DESeqPrep_fn(Counts)
TPM <- DESeqPrep_fn(TPM)
sizeFactors(Counts[[1]])
## NULL
sizeFactors(TPM[[1]])
## A1 A2 A3 B1 B2 B3
## 1.2059202 1.0600776 0.7425445 1.0155827 1.0482637 1.0759300
# Size factors don't exist in the Counts dds!
# Normalization factors are calculated in the Counts dds instead!
assays(Counts[[1]])
## List of length 6
## names(6): counts avgTxLength normalizationFactors mu H cooks
assays(TPM[[1]])
## List of length 4
## names(4): counts mu H cooks
colData(Counts[[1]])
## DataFrame with 6 rows and 3 columns
## Sample Group Path
## <factor> <factor> <character>
## A1 A1 A salmon_output/A1.sal..
## A2 A2 A salmon_output/A2.sal..
## A3 A3 A salmon_output/A3.sal..
## B1 B1 B salmon_output/B1.sal..
## B2 B2 B salmon_output/B2.sal..
## B3 B3 B salmon_output/B3.sal..
colData(TPM[[1]])
## DataFrame with 6 rows and 4 columns
## Sample Group Path sizeFactor
## <factor> <factor> <character> <numeric>
## A1 A1 A salmon_output/A1.sal.. 1.205920
## A2 A2 A salmon_output/A2.sal.. 1.060078
## A3 A3 A salmon_output/A3.sal.. 0.742545
## B1 B1 B salmon_output/B1.sal.. 1.015583
## B2 B2 B salmon_output/B2.sal.. 1.048264
## B3 B3 B salmon_output/B3.sal.. 1.075930
# Extract and save the size factors in a data frame
sizeFactor <- as.data.frame(round(sizeFactors(TPM[[1]]), 3))
colnames(sizeFactor) <- 'Size_Factor'
sizeFactor <- sizeFactor %>%
rownames_to_column(var="Sample") %>%
inner_join(metadata[, 1:ncol(metadata)-1], by="Sample")
# Create a plot comparing the size factors by sample
ggplot(sizeFactor, aes(x=Sample,
y=Size_Factor,
fill=Group,
label=Size_Factor)) +
geom_bar(stat="identity", width=0.8) +
theme_bw() +
ggtitle("Size Factors in TPM-DESeq") +
geom_text(vjust=1.5) +
theme(axis.text.x=element_text(angle=45,
vjust=0.5)) + ylab("Size Factor")
# Extract and normalization factors in a data frame
normf <- as.data.frame(normalizationFactors(Counts[[1]])) %>%
gather(Sample, Normalization_Factor) %>%
inner_join(metadata[, 1:2], by="Sample")
normf$Sample <- factor(normf$Sample, levels=SampleNames)
normf$Group <- factor(normf$Group, levels=GroupLevel)
# Create a scatter plot showing distribution of normalization factors
normFactors_plot <- ggplot(normf,
aes(x=Sample, y=Normalization_Factor)) +
geom_jitter(alpha=0.5, aes(color=Group)) +
# Add a boxplot to provide statistics in each sample
geom_boxplot(aes(x=Sample, y=Normalization_Factor),
outlier.shape=NA, alpha=0.5) +
theme_bw() +
ggtitle("Normalization Factors in Counts-DESeq") +
theme(axis.text.x=element_text(angle=45,
vjust=0.5)) +
ylab("Normalization Factor / Gene") +
# Add a dashed horizontal line to indicate where normalization factor=1
geom_hline(yintercept=1,
color="blue",
linetype="dashed")
# Print the normalization factor scatter plot
print(normFactors_plot)
# Print the same plot with larger y-magnification in order to observe the box plot
normFactors_plot +
ylim(0.5, 1.5)
# Assigne what to compare
GroupOfInterest <- Contrast[1]
# Set a function for a PCA plot
QCPCA_fn <- function(inputList, Title) {
plotPCA(inputList[[2]], # takes vsd
intgroup=GroupOfInterest) + theme_bw() + ggtitle(Title)
}
# Set heatmap annotation
ColOfInterest <- !colnames(metadata) %in% c("Sample", "Path")
HeatmapAnno <- as.data.frame(metadata[, ColOfInterest])
rownames(HeatmapAnno) <- SampleNames
colnames(HeatmapAnno) <- colnames(metadata)[ColOfInterest]
# Set a function for a correlation heatmap
QCcorrHeatmap_fn <- function(inputList, Title) {
# Extract transformed count matrix
mtx <- assay(inputList[[2]]) # takes vsd
# Calculate correlation and store in the matrix
mtx <- cor(mtx)
# Create a correlation heatmap
return(pheatmap(mtx,
annotation=HeatmapAnno,
main=paste("Sample Correlation Heatmap:",
Title)))
}
QCPCA_fn(TPM, "QC PCA: TPM")
QCPCA_fn(Counts, "QC PCA: Counts")
# TPM
QCcorrHeatmap_fn(TPM, "TPM")
# Counts
QCcorrHeatmap_fn(Counts, "Counts")
# Create a list for TPM and Counts dds
ddsList <- list(TPM=TPM[[1]], Counts=Counts[[1]])
for (x in TvC) {
# Run DESeq()
ddsList[[x]] <- DESeq(ddsList[[x]])
print(resultsNames(ddsList[[x]]))
}
## [1] "Intercept" "Group_B_vs_A"
## [1] "Intercept" "Group_B_vs_A"
# Plot dispersion
for (x in TvC) {
plotDispEsts(ddsList[[x]],
ylab="Dispersion",
xlab="Mean of Normalized Counts",
main=paste("Dispersion of", x, "Input"))
}
# Create an empty list for shrunken dds
shr_ddsList <- list(TPM=c(), Counts=c())
for (x in TvC) {
# shrink
shr_ddsList[[x]] <- lfcShrink(ddsList[[x]],
contrast=Contrast, # contrast
type="normal") # is paired with "normal" type
}
# Set FDR threshold
alpha=0.1
# FDR threshold vector
FDRv=c("< 0.1", "> 0.1")
# Set a function cleaning table
Sig_fn <- function(df, Input) {
df <- df %>%
rownames_to_column(var="Gene") %>%
mutate(FDR=ifelse(padj < 0.1 & !is.na(padj),
FDRv[1],
FDRv[2]),
Input=Input)
return(df)
}
# Initialize lists of result tables with (resList) or without (shr_resList) shrinkage
resList <- ddsList
shr_resList <- ddsList
for (x in TvC) {
# Extract results
resList[[x]] <- as.data.frame(results(ddsList[[x]],
contrast=Contrast,
alpha=alpha))
shr_resList[[x]] <- as.data.frame(shr_ddsList[[x]])
# clean the data frame
resList[[x]] <- Sig_fn(resList[[x]], x)
shr_resList[[x]] <- Sig_fn(shr_resList[[x]], x)
}
# No shrinkage summary
summary(resList)
## Length Class Mode
## TPM 9 data.frame list
## Counts 9 data.frame list
head(resList[['TPM']])
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 11.2413583 0.7644805 0.7163354 1.0672102 2.858769e-01
## 2 A3GALT2 3.7110877 4.2586719 1.6186278 2.6310384 8.512442e-03
## 3 A4GNT 3.8038346 1.2723597 1.1529779 1.1035421 2.697918e-01
## 4 AACSP1 55.7607695 -2.1138938 0.3648669 -5.7936029 6.889223e-09
## 5 AADACL2 0.1641094 0.8714154 4.0804729 0.2135575 8.308922e-01
## 6 AADACL4 0.5008675 -2.5075764 3.7818782 -0.6630506 5.072982e-01
## padj FDR Input
## 1 3.534571e-01 > 0.1 TPM
## 2 1.468366e-02 < 0.1 TPM
## 3 3.361484e-01 > 0.1 TPM
## 4 1.735453e-08 < 0.1 TPM
## 5 NA > 0.1 TPM
## 6 NA > 0.1 TPM
head(resList[['Counts']])
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 10.8066609 0.6941601 0.7188686 0.9656286 3.342301e-01
## 2 A3GALT2 3.8166660 4.2974513 1.6078205 2.6728426 7.521152e-03
## 3 A4GNT 3.8610666 1.2967288 1.1440446 1.1334600 2.570211e-01
## 4 AACSP1 56.5992344 -2.1461615 0.3638336 -5.8987446 3.662778e-09
## 5 AADACL2 0.1700555 0.9200020 4.0804729 0.2254645 8.216179e-01
## 6 AADACL4 0.5019710 -2.4653985 3.7632530 -0.6551243 5.123877e-01
## padj FDR Input
## 1 4.052284e-01 > 0.1 Counts
## 2 1.302526e-02 < 0.1 Counts
## 3 3.212300e-01 > 0.1 Counts
## 4 9.269316e-09 < 0.1 Counts
## 5 NA > 0.1 Counts
## 6 NA > 0.1 Counts
# Shrinkage summary
summary(shr_resList)
## Length Class Mode
## TPM 9 data.frame list
## Counts 9 data.frame list
head(shr_resList[['TPM']])
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 11.2413583 0.7473037 0.6999647 1.0672102 2.858769e-01
## 2 A3GALT2 3.7110877 3.7983247 1.4378564 2.6310384 8.512442e-03
## 3 A4GNT 3.8038346 1.2010948 1.0856049 1.1035421 2.697918e-01
## 4 AACSP1 55.7607695 -2.1013085 0.3625105 -5.7936029 6.889223e-09
## 5 AADACL2 0.1641094 0.4978588 2.3312198 0.2135575 8.308922e-01
## 6 AADACL4 0.5008675 -1.5560865 2.3037874 -0.6630506 5.072982e-01
## padj FDR Input
## 1 3.534571e-01 > 0.1 TPM
## 2 1.468366e-02 < 0.1 TPM
## 3 3.361484e-01 > 0.1 TPM
## 4 1.735453e-08 < 0.1 TPM
## 5 NA > 0.1 TPM
## 6 NA > 0.1 TPM
head(shr_resList[['Counts']])
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 10.8066609 0.6783316 0.7023183 0.9656286 3.342301e-01
## 2 A3GALT2 3.8166660 3.8371575 1.4296784 2.6728426 7.521152e-03
## 3 A4GNT 3.8610666 1.2249101 1.0778156 1.1334600 2.570211e-01
## 4 AACSP1 56.5992344 -2.1333964 0.3614677 -5.8987446 3.662778e-09
## 5 AADACL2 0.1700555 0.5246799 2.3270647 0.2254645 8.216179e-01
## 6 AADACL4 0.5019710 -1.5324064 2.2975661 -0.6551243 5.123877e-01
## padj FDR Input
## 1 4.052284e-01 > 0.1 Counts
## 2 1.302526e-02 < 0.1 Counts
## 3 3.212300e-01 > 0.1 Counts
## 4 9.269316e-09 < 0.1 Counts
## 5 NA > 0.1 Counts
## 6 NA > 0.1 Counts
# Set ylim: has to adjusted by users depending on data
yl <- c(-20, 20)
# Set min log2 fold change of interest
mLog <- c(-1, 1)
# Define a function creating an MA plot
MA_fn <- function(List, Shr) {
MAList <- ddsList
for (i in 1:2) {
MAplot <- ggplot(List[[i]],
aes(x=baseMean,
y=log2FoldChange,
color=FDR)) + geom_point() + scale_x_log10() + theme_bw() + scale_color_manual(values=c("blue", "grey")) + ggtitle(paste("MA plot:", names(List)[i], "Input with", Shr)) + ylim(yl[1], yl[2])+ geom_hline(yintercept=c(mLog[1], mLog[2]), linetype="dashed", color="red")
MAList[[i]] <- MAplot
}
return(MAList)
}
# Create MA plots with or without shrinkage and store in a list
MA <- MA_fn(resList, "No Shrinkage")
shr_MA <- MA_fn(shr_resList, "Shrinkage")
# TPM with or without shrinkage
grid.arrange(MA[[1]], shr_MA[[1]], nrow=1)
# TPM with or without shrinkage
grid.arrange(MA[[2]], shr_MA[[2]], nrow=1)
# Combining total data table
res <- rbind(resList[['TPM']], resList[['Counts']])
res$Input <- factor(res$Input, levels=TvC) # TvC=c("TPM", "Counts")
# Create a plot presenting distribution of FDR
ggplot(res,
aes(x=padj, color=Input)) +
geom_density(size=1, aes(y=..count..)) +
theme_bw() +
ggtitle("Distribution of False Discovery Rate (FDR)") +
xlab("Adjusted P-Value (FDR)") +
ylab("Count") +
geom_vline(xintercept=alpha,
color="black",
linetype="dashed",
size=1) +
scale_x_continuous(breaks=seq(0, 1, by=0.1))
ggplot(res[res$FDR == "< 0.1", ], # Subset rows with FDR < alpha
aes(x=log2FoldChange,
color=Input)) +
geom_density(size=1, aes(y=..count..)) +
theme_bw() +
ylab("Count") +
geom_vline(xintercept=c(mLog[1], mLog[2]),
color="black",
linetype="dashed", size=1) +
ggtitle("Distribution of Log2 Folds by Input Type") +
xlim(xlim[1], xlim[2])
# Initialize a list
lowfdrList <- ddsList # A list for normalized counts matrix with FDR below alpha
highfoldList <- ddsList # A list for normalized counts with log2foldchange over minLog
for (x in TvC) {
# Create filtering vectors with alpha and log2foldchange
BelowAlpha <- resList[[x]]$FDR == FDRv[1]
overmLog <- abs(resList[[x]]$log2FoldChange) > mLog[2] # mLog has been set to c(-1, 1) previously
# Extract transformed counts from vsd objects (TPM[['vsd']] and Counts[['vsd']])
if (x == "TPM") {
normCounts <- counts(TPM[['dds']], normalized=T)
} else {
normCounts <- counts(Counts[['dds']], normalized=T)
}
# Update the normalized count matrix with FDR below alpha
lowfdrList[[x]] <- normCounts[BelowAlpha, ]
highfoldList[[x]] <- normCounts[BelowAlpha & overmLog, ]
summary(lowfdrList[[x]])
summary(highfoldList[[x]])
}
# Initialize map lists
lowfdrMap <- ddsList
highfoldMap <- ddsList
# Set a function creating a heatmap
ProfileHeatmap_fn <- function(inputmatrix, Title1, Title2, Title3=NULL) {
as.ggplot(pheatmap(inputmatrix,
annotation=HeatmapAnno,
scale="row", # presents z-score instead of counts
show_rownames=F,
main=paste("Transcription Profiles with",
Title1,
"input and",
Title2,
alpha,
Title3)))
}
# Create and save heatmaps
for (x in TvC) {
lowfdrMap[[x]] <- ProfileHeatmap_fn(lowfdrList[[x]],
Title1=x,
Title2="FDR <")
highfoldMap[[x]] <- ProfileHeatmap_fn(highfoldList[[x]],
Title1=x,
Title2="FDR <",
Title3=paste("+ Absolte Log2 Fold Change >", mLog[2]))
}
# Count number of NA genes
type=c("Zero Counts", "Outliers", "Total NA Genes")
NAstat <- res %>%
group_by(Input) %>%
summarize(zero=sum(is.na(log2FoldChange)),
outlier=sum(is.na(pvalue) & is.na(padj))) %>%
mutate(total=sum(zero, outlier)) %>%
gather(Type, Number, -Input) %>%
mutate(Type=factor(case_when(Type == "zero" ~ type[1],
Type == "outlier" ~ type[2],
Type == "total" ~ type[3]),
levels=type))
# Plot number of NA genes
ggplot(NAstat, aes(x=Type, y=Number, group=Input, fill=Input, label=Number)) +
geom_bar(stat="identity", position="dodge") +
theme_bw() +
geom_text(position=position_dodge(width=1), vjust=1.5) +
ggtitle("Number of NA Genes") +
ylab("Number of Genes")
FDRrankList: ranked by FDR
lfcList: ranked by absolute fold change
UPlfcList: ranked by magnitude of fold change increase
DOWNlfcList: ranked by manitude of fold change decrease
# Create a new list having DE table with FDR below alpha
lowfdr_resList <- resList
for (x in TvC) {
lowfdr_resList[[x]] <- filter(resList[[x]],
FDR == FDRv[1]) %>%
as.data.table()
}
# Initialize new lists in order to store rank-updated result DE tables
FDRrankList <- lowfdr_resList
lfcList <- lowfdr_resList
UPlfcList <- lowfdr_resList
DOWNlfcList <- lowfdr_resList
# Set a function creating a column for the rank
Ranking_fn <- function(x) {mutate(x, Rank=1:nrow(x))}
for (x in TvC) {
# Rearrange genes with FDR
FDRrankList[[x]] <- lowfdr_resList[[x]][order(padj),] %>%
Ranking_fn()
# Rearrange genew with absolute log2FoldChange
lfcList[[x]] <- lowfdr_resList[[x]][order(-abs(log2FoldChange)),] %>%
Ranking_fn()
# Rearrange genes with log2FoldChange (decreasing order)
UPlfcList[[x]] <- lowfdr_resList[[x]][order(-log2FoldChange),] %>%
Ranking_fn()
# Rearrange genes with log2FoldChange (increasing order)
DOWNlfcList[[x]] <- lowfdr_resList[[x]][order(log2FoldChange),] %>%
Ranking_fn()
}
# Explore the ranks
print(c(FDRrankList, lfcList, UPlfcList, DOWNlfcList))
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: ACBD7 3168.682742 2.9056054 0.06424347 45.228029 0.00000000
## 2: AFP 32584.595366 -11.0421341 0.15395409 -71.723551 0.00000000
## 3: AKAP11 13905.676097 2.4997058 0.03666677 68.173600 0.00000000
## 4: AMIGO1 3568.001961 4.4114268 0.06849514 64.404960 0.00000000
## 5: AMY2B 1557.113021 3.6430792 0.08662900 42.053805 0.00000000
## ---
## 4726: SALL4P1 1.454384 -4.0531916 2.21570166 -1.829304 0.06735410
## 4727: GDF9 70.162940 0.5210032 0.28544731 1.825217 0.06796832
## 4728: RPL12P10 3.588481 -2.6723591 1.46472861 -1.824474 0.06808046
## 4729: RPL34P18 9.116960 -1.4682061 0.80491820 -1.824044 0.06814545
## 4730: IMPDH1P10 12.307808 -7.1302067 3.91051047 -1.823344 0.06825127
## padj FDR Input Rank
## 1: 0.00000000 < 0.1 TPM 1
## 2: 0.00000000 < 0.1 TPM 2
## 3: 0.00000000 < 0.1 TPM 3
## 4: 0.00000000 < 0.1 TPM 4
## 5: 0.00000000 < 0.1 TPM 5
## ---
## 4726: 0.09876745 < 0.1 TPM 4726
## 4727: 0.09964469 < 0.1 TPM 4727
## 4728: 0.09978798 < 0.1 TPM 4728
## 4729: 0.09986212 < 0.1 TPM 4729
## 4730: 0.09999605 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: ABCD2 1069.176400 5.560740 0.14812494 37.540872 0.00000000
## 2: ACBD7 3223.871084 2.902042 0.06325036 45.881828 0.00000000
## 3: AFP 32279.951238 -11.043880 0.14081349 -78.429135 0.00000000
## 4: AKAP11 14142.443732 2.494980 0.03650368 68.348719 0.00000000
## 5: AMIGO1 3630.556248 4.406517 0.06828149 64.534572 0.00000000
## ---
## 4730: SUGT1P2 10.591190 -1.285347 0.70388149 -1.826085 0.06783752
## 4731: MTND2P5 4.330142 2.185290 1.19716309 1.825390 0.06794219
## 4732: RPSAP21 2.057284 -3.437381 1.88421063 -1.824308 0.06810546
## 4733: OR52W1 3.905048 2.282493 1.25138546 1.823973 0.06815620
## 4734: SALL4P1 1.468341 -4.021370 2.20511483 -1.823656 0.06820416
## padj FDR Input Rank
## 1: 0.00000000 < 0.1 Counts 1
## 2: 0.00000000 < 0.1 Counts 2
## 3: 0.00000000 < 0.1 Counts 3
## 4: 0.00000000 < 0.1 Counts 4
## 5: 0.00000000 < 0.1 Counts 5
## ---
## 4730: 0.09937551 < 0.1 Counts 4730
## 4731: 0.09950780 < 0.1 Counts 4731
## 4732: 0.09972586 < 0.1 Counts 4732
## 4733: 0.09977907 < 0.1 Counts 4733
## 4734: 0.09982819 < 0.1 Counts 4734
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1935.1129 -14.42718266 1.18067729 -12.219412 2.448482e-34
## 2: GAGE12C 665.4894 -12.88787310 1.18254611 -10.898411 1.172881e-27
## 3: GAGE13 612.2774 -12.76494282 1.18202949 -10.799175 3.473116e-27
## 4: GAGE12F 2805.4936 -12.51900730 0.83651932 -14.965593 1.232215e-50
## 5: UTP14C 531.7746 12.47219039 1.18872293 10.492092 9.392533e-26
## ---
## 4726: GM2A 5597.9316 0.10758757 0.04117832 2.612724 8.982382e-03
## 4727: CD2BP2 6723.4716 -0.10030164 0.04900183 -2.046896 4.066833e-02
## 4728: UQCRFS1 5731.6147 -0.09801823 0.04466793 -2.194376 2.820840e-02
## 4729: PRPF6 8832.8778 -0.09437979 0.04835202 -1.951931 5.094644e-02
## 4730: SELENOP 20675.0072 -0.08968953 0.04428648 -2.025212 4.284557e-02
## padj FDR Input Rank
## 1: 1.092594e-33 < 0.1 TPM 1
## 2: 4.687467e-27 < 0.1 TPM 2
## 3: 1.378505e-26 < 0.1 TPM 3
## 4: 6.681727e-50 < 0.1 TPM 4
## 5: 3.642432e-25 < 0.1 TPM 5
## ---
## 4726: 1.543081e-02 < 0.1 TPM 4726
## 4727: 6.294319e-02 < 0.1 TPM 4727
## 4728: 4.503207e-02 < 0.1 TPM 4728
## 4729: 7.715446e-02 < 0.1 TPM 4729
## 4730: 6.601152e-02 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1975.0482 -14.41206728 1.18066330 -12.206755 2.860740e-34
## 2: GAGE12C 679.3481 -12.87323165 1.18248182 -10.886621 1.335017e-27
## 3: GAGE13 624.8986 -12.74994147 1.18192899 -10.787401 3.947987e-27
## 4: GAGE12F 2863.4359 -12.51268995 0.83641254 -14.959950 1.341265e-50
## 5: UTP14C 541.2130 12.50175251 1.18864170 10.517680 7.161562e-26
## ---
## 4730: CD2BP2 6800.0952 -0.10470103 0.04950383 -2.115009 3.442919e-02
## 4731: GM2A 5697.3428 0.10317346 0.04064652 2.538310 1.113893e-02
## 4732: UQCRFS1 5839.0452 -0.10252336 0.04427907 -2.315391 2.059153e-02
## 4733: PRPF6 8994.9911 -0.09892368 0.04794184 -2.063410 3.907370e-02
## 4734: SELENOP 21049.9178 -0.09447398 0.04403219 -2.145566 3.190761e-02
## padj FDR Input Rank
## 1: 1.266586e-33 < 0.1 Counts 1
## 2: 5.304090e-27 < 0.1 Counts 2
## 3: 1.554295e-26 < 0.1 Counts 3
## 4: 7.238026e-50 < 0.1 Counts 4
## 5: 2.758336e-25 < 0.1 Counts 5
## ---
## 4730: 5.388748e-02 < 0.1 Counts 4730
## 4731: 1.887083e-02 < 0.1 Counts 4731
## 4732: 3.361892e-02 < 0.1 Counts 4732
## 4733: 6.051445e-02 < 0.1 Counts 4733
## 4734: 5.024724e-02 < 0.1 Counts 4734
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: UTP14C 531.7746 12.47219 1.1887229 10.49209 9.392533e-26
## 2: NEUROD6 914.0304 12.29159 1.1806402 10.41096 2.209892e-25
## 3: CARTPT 444.8165 12.21474 1.1829535 10.32563 5.396390e-25
## 4: HBB 2497.0069 11.87435 0.7237799 16.40603 1.731676e-60
## 5: SLC17A6 641.0675 11.78002 1.1827184 9.96012 2.277885e-23
## ---
## 4726: BPIFB2 1748.3190 -12.42289 1.0233841 -12.13902 6.560335e-34
## 4727: GAGE12F 2805.4936 -12.51901 0.8365193 -14.96559 1.232215e-50
## 4728: GAGE13 612.2774 -12.76494 1.1820295 -10.79917 3.473116e-27
## 4729: GAGE12C 665.4894 -12.88787 1.1825461 -10.89841 1.172881e-27
## 4730: RHOXF2B 1935.1129 -14.42718 1.1806773 -12.21941 2.448482e-34
## padj FDR Input Rank
## 1: 3.642432e-25 < 0.1 TPM 1
## 2: 8.512813e-25 < 0.1 TPM 2
## 3: 2.060440e-24 < 0.1 TPM 3
## 4: 1.031859e-59 < 0.1 TPM 4
## 5: 8.459668e-23 < 0.1 TPM 5
## ---
## 4726: 2.904992e-33 < 0.1 TPM 4726
## 4727: 6.681727e-50 < 0.1 TPM 4727
## 4728: 1.378505e-26 < 0.1 TPM 4728
## 4729: 4.687467e-27 < 0.1 TPM 4729
## 4730: 1.092594e-33 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: UTP14C 541.2130 12.50175 1.1886417 10.517680 7.161562e-26
## 2: NEUROD6 929.8364 12.30527 1.1805725 10.423137 1.944330e-25
## 3: CARTPT 452.4859 12.24355 1.1828268 10.351089 4.137522e-25
## 4: HBB 2400.0112 12.13599 0.8113413 14.957939 1.382416e-50
## 5: SLC17A6 652.2483 11.81567 1.1825758 9.991472 1.660961e-23
## ---
## 4730: BPIFB2 1782.0106 -12.45487 1.0233020 -12.171251 4.422530e-34
## 4731: GAGE12F 2863.4359 -12.51269 0.8364125 -14.959950 1.341265e-50
## 4732: GAGE13 624.8986 -12.74994 1.1819290 -10.787401 3.947987e-27
## 4733: GAGE12C 679.3481 -12.87323 1.1824818 -10.886621 1.335017e-27
## 4734: RHOXF2B 1975.0482 -14.41207 1.1806633 -12.206755 2.860740e-34
## padj FDR Input Rank
## 1: 2.758336e-25 < 0.1 Counts 1
## 2: 7.430922e-25 < 0.1 Counts 2
## 3: 1.574349e-24 < 0.1 Counts 3
## 4: 7.454287e-50 < 0.1 Counts 4
## 5: 6.131486e-23 < 0.1 Counts 5
## ---
## 4730: 1.953073e-33 < 0.1 Counts 4730
## 4731: 7.238026e-50 < 0.1 Counts 4731
## 4732: 1.554295e-26 < 0.1 Counts 4732
## 4733: 5.304090e-27 < 0.1 Counts 4733
## 4734: 1.266586e-33 < 0.1 Counts 4734
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1935.1129 -14.42718 1.1806773 -12.21941 2.448482e-34
## 2: GAGE12C 665.4894 -12.88787 1.1825461 -10.89841 1.172881e-27
## 3: GAGE13 612.2774 -12.76494 1.1820295 -10.79917 3.473116e-27
## 4: GAGE12F 2805.4936 -12.51901 0.8365193 -14.96559 1.232215e-50
## 5: BPIFB2 1748.3190 -12.42289 1.0233841 -12.13902 6.560335e-34
## ---
## 4726: SLC17A6 641.0675 11.78002 1.1827184 9.96012 2.277885e-23
## 4727: HBB 2497.0069 11.87435 0.7237799 16.40603 1.731676e-60
## 4728: CARTPT 444.8165 12.21474 1.1829535 10.32563 5.396390e-25
## 4729: NEUROD6 914.0304 12.29159 1.1806402 10.41096 2.209892e-25
## 4730: UTP14C 531.7746 12.47219 1.1887229 10.49209 9.392533e-26
## padj FDR Input Rank
## 1: 1.092594e-33 < 0.1 TPM 1
## 2: 4.687467e-27 < 0.1 TPM 2
## 3: 1.378505e-26 < 0.1 TPM 3
## 4: 6.681727e-50 < 0.1 TPM 4
## 5: 2.904992e-33 < 0.1 TPM 5
## ---
## 4726: 8.459668e-23 < 0.1 TPM 4726
## 4727: 1.031859e-59 < 0.1 TPM 4727
## 4728: 2.060440e-24 < 0.1 TPM 4728
## 4729: 8.512813e-25 < 0.1 TPM 4729
## 4730: 3.642432e-25 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1975.0482 -14.41207 1.1806633 -12.206755 2.860740e-34
## 2: GAGE12C 679.3481 -12.87323 1.1824818 -10.886621 1.335017e-27
## 3: GAGE13 624.8986 -12.74994 1.1819290 -10.787401 3.947987e-27
## 4: GAGE12F 2863.4359 -12.51269 0.8364125 -14.959950 1.341265e-50
## 5: BPIFB2 1782.0106 -12.45487 1.0233020 -12.171251 4.422530e-34
## ---
## 4730: SLC17A6 652.2483 11.81567 1.1825758 9.991472 1.660961e-23
## 4731: HBB 2400.0112 12.13599 0.8113413 14.957939 1.382416e-50
## 4732: CARTPT 452.4859 12.24355 1.1828268 10.351089 4.137522e-25
## 4733: NEUROD6 929.8364 12.30527 1.1805725 10.423137 1.944330e-25
## 4734: UTP14C 541.2130 12.50175 1.1886417 10.517680 7.161562e-26
## padj FDR Input Rank
## 1: 1.266586e-33 < 0.1 Counts 1
## 2: 5.304090e-27 < 0.1 Counts 2
## 3: 1.554295e-26 < 0.1 Counts 3
## 4: 7.238026e-50 < 0.1 Counts 4
## 5: 1.953073e-33 < 0.1 Counts 5
## ---
## 4730: 6.131486e-23 < 0.1 Counts 4730
## 4731: 7.454287e-50 < 0.1 Counts 4731
## 4732: 1.574349e-24 < 0.1 Counts 4732
## 4733: 7.430922e-25 < 0.1 Counts 4733
## 4734: 2.758336e-25 < 0.1 Counts 4734
# Set a function rebuilding DE tables with gene ranks
combineTable_fn <- function(List){
# Select columns and join the data frames by gene
full_join(List[[TvC[1]]][,.(Gene, Input, Rank, baseMean)],
List[[TvC[2]]][,.(Gene, Input, Rank, baseMean)], by="Gene") %>%
# Add columns assining gene expression levels, rank differences, and mean ranks
mutate(logMeanExpression=log(baseMean.x+baseMean.y/2),
RankDiff=Rank.x-Rank.y,
MeanRank=(Rank.x+Rank.y)/2)
}
# Explore outputs of the function
head(combineTable_fn(FDRrankList))
## Gene Input.x Rank.x baseMean.x Input.y Rank.y baseMean.y logMeanExpression
## 1: ACBD7 TPM 1 3168.683 Counts 2 3223.871 8.472325
## 2: AFP TPM 2 32584.595 Counts 3 32279.951 10.793939
## 3: AKAP11 TPM 3 13905.676 Counts 4 14142.444 9.951177
## 4: AMIGO1 TPM 4 3568.002 Counts 5 3630.556 8.591053
## 5: AMY2B TPM 5 1557.113 Counts 6 1581.779 7.761320
## 6: ANP32B TPM 6 10924.599 Counts 7 11143.790 9.710903
## RankDiff MeanRank
## 1: -1 1.5
## 2: -1 2.5
## 3: -1 3.5
## 4: -1 4.5
## 5: -1 5.5
## 6: -1 6.5
dim(combineTable_fn(FDRrankList))
## [1] 4761 10
tail(combineTable_fn(FDRrankList))
## Gene Input.x Rank.x baseMean.x Input.y Rank.y baseMean.y
## 1: INSL6 <NA> NA NA Counts 4723 8.692921
## 2: ASS1P10 <NA> NA NA Counts 4726 9.072376
## 3: RPS3AP30 <NA> NA NA Counts 4729 1.481576
## 4: SUGT1P2 <NA> NA NA Counts 4730 10.591190
## 5: MTND2P5 <NA> NA NA Counts 4731 4.330142
## 6: OR52W1 <NA> NA NA Counts 4733 3.905048
## logMeanExpression RankDiff MeanRank
## 1: NA NA NA
## 2: NA NA NA
## 3: NA NA NA
## 4: NA NA NA
## 5: NA NA NA
## 6: NA NA NA
# Set a function plotting gene ranks between TPM- (x-axis) and Counts-Inputs (y-axis)
compareRanks_fn <- function(df, rankedby) {
ggplot(df,
aes(x=Rank.x,
y=Rank.y,
color=logMeanExpression)) +
geom_point(alpha=0.5) +
theme_bw() +
xlab("Rank with TPM") +
ylab("Rank with Counts") +
geom_abline(slope=1, color="black", size=0.5) +
ggtitle(paste(rankedby, "Ranking with TPM vs Counts Inputs")) +
scale_color_gradient(low="blue", high="red")
}
# Set a function plotting the rank difference over the gene expression level
RankdiffOverMean_fn <- function(df, rankedby) {
ggplot(df,
aes(x=logMeanExpression,
y=RankDiff,
color=MeanRank)) +
geom_point(alpha=0.5) +
theme_bw() +
ylab("Rank Difference (TPM - Counts)") +
ggtitle(paste("Rank Difference Inputs (TPM - Counts) in", rankedby)) +
geom_hline(yintercept=0, color="black", size=0.5) +
scale_color_gradient(low="blue", high="red")
}
# Ranked by FDR
compareRanks_fn(combineTable_fn(FDRrankList),
"FDR")
# Ranked by absolute fold change
compareRanks_fn(combineTable_fn(lfcList),
"Absolute Log2FoldChange")
# Ranked by magnitude of positive fold change
compareRanks_fn(combineTable_fn(UPlfcList),
"Log2FoldChange (Increased)")
# Ranked by magnitude of negative fold change
compareRanks_fn(combineTable_fn(DOWNlfcList),
"Log2FoldChange (Decreased)")
# Ranked by FDR
RankdiffOverMean_fn(combineTable_fn(FDRrankList),
"FDR")
# Ranked by absolute fold change
RankdiffOverMean_fn(combineTable_fn(lfcList),
"Absolute Log2FoldChange")
# Ranked by magnitude of positive fold change
RankdiffOverMean_fn(combineTable_fn(UPlfcList),
"Log2FoldChange (Increased)")
# Ranked by magnitude of negative fold change
RankdiffOverMean_fn(combineTable_fn(DOWNlfcList),
"Log2FoldChange (Decreased)")
# Clean data to generate an upset plot
res <- res %>%
# Filter genes with valid padj
filter(!is.na(padj)) %>%
# Detect genes which are up/down/unchanged change patterns in either TPM and Counts inputs
mutate(Up=ifelse(FDR == FDRv[1] & log2FoldChange > 0, Gene, ""), # What are upregulated genes?
Down=ifelse(FDR == FDRv[1] & log2FoldChange < 0, Gene, ""), # What are downregulated genes?
Unchanged=ifelse(FDR == FDRv[2], Gene, ""), # What are unchanged genes?
TPM_Input=ifelse(Input == "TPM", Gene, ""), # What are the genes from TPM input?
Counts_Input=ifelse(Input == "Counts", Gene, "")) # What are the genes from Counts input?
# Create a list storing groups of interest
upsetInput <- list(Up=res$Up,
Down=res$Down,
Unchanged=res$Unchanged,
TPM_Input=res$TPM,
Counts_Input=res$Counts)
# Create the upset plot
upset(fromList(upsetInput),
sets.x.label="Number of Genes per Group",
order.by="freq")
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/mira/miniconda3/envs/r/lib/libopenblasp-r0.3.10.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] AnnotationDbi_1.52.0 UpSetR_1.4.0
## [3] ggrepel_0.8.2 ggplotify_0.0.5
## [5] gridExtra_2.3 pheatmap_1.0.12
## [7] DESeq2_1.30.0 SummarizedExperiment_1.20.0
## [9] Biobase_2.50.0 MatrixGenerics_1.2.0
## [11] matrixStats_0.57.0 GenomicRanges_1.42.0
## [13] GenomeInfoDb_1.26.0 IRanges_2.24.0
## [15] S4Vectors_0.28.0 tximport_1.18.0
## [17] forcats_0.5.0 stringr_1.4.0
## [19] dplyr_1.0.2 purrr_0.3.4
## [21] readr_1.4.0 tidyr_1.1.2
## [23] tibble_3.0.4 ggplot2_3.3.2
## [25] tidyverse_1.3.0 AnnotationHub_2.22.0
## [27] BiocFileCache_1.14.0 dbplyr_2.0.0
## [29] BiocGenerics_0.36.0 rmarkdown_2.5
## [31] data.table_1.13.2
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1
## [3] XVector_0.30.0 fs_1.5.0
## [5] rstudioapi_0.13 farver_2.0.3
## [7] bit64_4.0.5 interactiveDisplayBase_1.28.0
## [9] fansi_0.4.1 lubridate_1.7.9.2
## [11] xml2_1.3.2 splines_4.0.3
## [13] geneplotter_1.68.0 knitr_1.30
## [15] jsonlite_1.7.1 broom_0.7.2
## [17] annotate_1.68.0 shiny_1.5.0
## [19] BiocManager_1.30.10 compiler_4.0.3
## [21] httr_1.4.2 rvcheck_0.1.8
## [23] backports_1.2.0 assertthat_0.2.1
## [25] Matrix_1.2-18 fastmap_1.0.1
## [27] cli_2.2.0 later_1.1.0.1
## [29] htmltools_0.5.0 tools_4.0.3
## [31] gtable_0.3.0 glue_1.4.2
## [33] GenomeInfoDbData_1.2.4 rappdirs_0.3.1
## [35] Rcpp_1.0.5 cellranger_1.1.0
## [37] vctrs_0.3.5 xfun_0.19
## [39] ps_1.5.0 rvest_0.3.6
## [41] mime_0.9 lifecycle_0.2.0
## [43] XML_3.99-0.3 zlibbioc_1.36.0
## [45] scales_1.1.1 hms_0.5.3
## [47] promises_1.1.1 RColorBrewer_1.1-2
## [49] yaml_2.2.1 curl_4.3
## [51] memoise_1.1.0 stringi_1.5.3
## [53] RSQLite_2.2.1 BiocVersion_3.12.0
## [55] genefilter_1.72.0 BiocParallel_1.24.0
## [57] rlang_0.4.9 pkgconfig_2.0.3
## [59] bitops_1.0-6 evaluate_0.14
## [61] lattice_0.20-41 labeling_0.4.2
## [63] bit_4.0.4 tidyselect_1.1.0
## [65] plyr_1.8.6 magrittr_2.0.1
## [67] R6_2.5.0 generics_0.1.0
## [69] DelayedArray_0.16.0 DBI_1.1.0
## [71] pillar_1.4.7 haven_2.3.1
## [73] withr_2.3.0 survival_3.2-7
## [75] RCurl_1.98-1.2 modelr_0.1.8
## [77] crayon_1.3.4 utf8_1.1.4
## [79] locfit_1.5-9.4 grid_4.0.3
## [81] readxl_1.3.1 blob_1.2.1
## [83] reprex_0.3.0 digest_0.6.27
## [85] xtable_1.8-4 httpuv_1.5.4
## [87] gridGraphics_0.5-0 munsell_0.5.0